clear;close all;
load data.mat;

seismap = interp1([0, 1, 32, 64], [1,1,1; 0,0,0;1,1,1;1,0,0], linspace(-1, 1, 65).^3*32+32);
%% fold number
mp = arrayfun(@(x) (x.SourceX + x.GroupX)/2, th);
mplist = unique(sort(mp));
foldnum = zeros(size(mplist));
for i = 1:length(mplist)
    foldnum(i) = sum(mp==mplist(i));
end
figure;
stem(mplist, foldnum);
%% common middle gather
[data, trhdr] = getcmg(dat, th, 20000);
plotseiswave(data, [trhdr.offset], dt, 'gather', 'line');

%% common shot gather
[data, trhdr] = getcsg(dat, th, 200*30);
plotseiswave(data, [trhdr.offset], dt, 'gather');
%%
dat = detrend(dat, 'linear');
%% balance
% dat = dat./repmat(max(abs(dat), [], 1), size(dat, 1), 1);
%% fk
close all;
% sourceX = [th.SourceX];
% shotX = unique(sourceX);
% for shotnum = shotX
%     trcID = find(sourceX==shotnum);
%     fdata = fkfilter(dat(:, trcID), th(trcID), dt, 3700);
%     dat(:, trcID) = fdata;
% end
% [data, trhdr] = getcsg(dat, th, 200*20);
% fdata = fkfilter(data, trhdr, dt, 3700);
% plotseiswave(data, [trhdr.offset], dt, 'trace');colormap(seismap);caxis([-1, 1]*0.5);
% plotseiswave(fdata, [trhdr.offset], dt, 'trace');colormap(seismap);caxis([-1, 1]*0.5);

%% compensation
% [data, trhdr] = getcmg(dat, th, 20000);
% plotseiswave(data, [trhdr.offset], dt, 'trace', 'line');
% gdata = gain(data, trhdr, dt, 0.0, 1.01, 'power');
% plotseiswave(gdata, [trhdr.offset], dt, 'trace', 'line');
% dat = gain(dat, th, dt, 0.0, 1.01, 'power');
%% SD deconvolution
close all;
[data, trhdr] = getcsg(dat, th, 200*20);
plotseiswave(data, [trhdr.offset], dt, 'trace', 'image');colormap(seismap);caxis([-1, 1]*0.5);
ds = spike_decon(data, 0.2, 100, dt);
plotseiswave(ds, [trhdr.offset], dt, 'trace', 'image');colormap(seismap);caxis([-1, 1]*0.5);
ds2 = zeros(size(data));
for i = 1:size(data, 2)
    ds2(:, i) = spike_decon(data(:, i), 0.4, 60, dt);
end
plotseiswave(ds2, [trhdr.offset], dt, 'trace', 'image');colormap(seismap);caxis([-1, 1]*0.5);
[b, a] = butter(4, [5, 20]*2*dt);
ds3 = filtfilt(b, a, ds);
plotseiswave(ds3, [trhdr.offset], dt, 'trace', 'image');colormap(seismap);caxis([-1, 1]*0.5);
plotseiswave(ds3, [trhdr.offset], dt, 'trace', 'line');

%% predictive deconv
close all;
[data, trhdr] = getcmg(dat, th, 200*20);
plotseiswave(data, [trhdr.offset], dt, 'trace', 'line');colormap(seismap);caxis([-1, 1]*0.5);
ds = zeros(size(data));
for i = 1:size(data, 2)
    [~, ds(:, i)] = predictive(data(:, i), size(data, 1), round(0.5/dt), 0.1);
end
plotseiswave(ds, [trhdr.offset], dt, 'trace', 'line');colormap(seismap);caxis([-1, 1]*0.5);

%% velocity analyse
Npicks = 100;
idxs = round((0:Npicks+1)/(Npicks+1)*(length(mplist)-1))+1;
idxs = idxs(2:Npicks+1);
vels = zeros(size(dat, 1), Npicks);
middlePoint = mplist(idxs);
for i = 74:Npicks
    close all;
    idx = idxs(i);
    fprintf('progress: %d/%d\n', i, Npicks);
    [data, trhdr] = getcmg(dat, th, mplist(idx));
    [vat, vav] = velocityanalyse(data, trhdr, dt, [1000, 6000], 50, 50, 0.1, 1);
    vels(:, i) = fillvel(vat, vav, (0:size(data, 1)-1)*dt);
end
close all;
%%
% [data, trhdr] = getcmg(dat, th, 20000);
% VelocityAnalysis(data, (0:size(data, 1)-1)*dt, [trhdr.offset], 500, 6000, (6000-500)/20,0,0,2,[3,3])